
/***************************************************************
 * sep_augxyz.c
 * Merge two xyz files: Particles given in the second xyz-file 
 * is stacked on top of the particles given the first. 
 * Writes new xyz file in start.xyz 
 * 
 * usage: sep_augxyz <xyz-file 1> <xyz-file 2> <dist (optional)> 
 *
 * Author: Jesper Hansen
 ***************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define VERSION "0.1"


void printpre( void);
void merror(const char *str, const char *fun, int line);
void readxyzfile(double *lbox_file, int *npart_file, const char *file);
void writexyzfile(double *lbox_1, double *lbox_2, int npart_1, int npart_2, 
		  char *file_1 , char *file_2, double dz);
void checkfile( void );


int main(int argc, char **argv){
  double lbox_1[3], lbox_2[3];
  int npart_1, npart_2;
 
  if ( argc < 3 || argc > 4 ){
    fprintf(stderr, 
	    "Usage is:  sep_augxyz <xyz-file 1> <xyz-file 2> <dist (optional)> \n");
    exit(EXIT_FAILURE);
  }

  if ( strcmp(argv[2], "start.xyz") == 0 ){
    fprintf(stderr, "Warning: ");
    fprintf(stderr, "Input argument is start.xyz wierd things will happen!\n");
  }

  FILE *fout = fopen("start.xyz", "r");
  if ( fout != NULL ){
    fprintf(stderr, "Warning: ");
    fprintf(stderr, "start.xyz exists - it will be overwritten\n");
    fclose(fout);
  }

  printpre();
  checkfile();

  double dz = 1.0;
  if ( argc==4 ) dz = atof(argv[3]);

  readxyzfile(lbox_1, &npart_1, argv[1]);
  readxyzfile(lbox_2, &npart_2, argv[2]);

  writexyzfile(lbox_1, lbox_2, npart_1, npart_2, argv[1], argv[2], dz);

  
  return 0;

}


void checkfile( void ){

  FILE *fout = fopen("start.xyz", "r");
  if ( fout != NULL ){
    fprintf(stderr, "Warning: ");
    fprintf(stderr, "start.xyz exists - it will be overwritten\n");
    fclose(fout);
  }

}



void printpre( void){

  printf("This is sep_augxyz version %s\n", VERSION);
  printf("This program is a part of the sep library - there is absolutely ");
  printf("NO WARRANTY\n\n");

}

void merror(const char *str, const char *fun, int line){

  fprintf(stderr, "%s at line %d, %s BAILING OUT\n", fun, line, str);

  exit(EXIT_FAILURE);

}




void readxyzfile(double *lbox_file, int *npart_file, const char *file){
  
  FILE *fin = fopen(file, "r");
  if ( fin==NULL )
    merror("%s at %d: Couldn't open file", __func__, __LINE__);
  
  if ( fscanf(fin, "%d\n", npart_file) != 1 )
    merror("%s at %d: Reading failure", __func__, __LINE__);

  if ( fscanf(fin, "%lf%lf%lf\n", &(lbox_file[0]), 
	      &(lbox_file[1]), &(lbox_file[2])) != 3 )
    merror("%s at %d: Reading failure", __func__, __LINE__);
       
  fclose(fin);

  printf("Read file %s\n", file);
  printf("Found %d particles and read box lengths %.3f %.3f %.3f\n", 
	 *npart_file, lbox_file[0], lbox_file[1], lbox_file[2]);

}



void writexyzfile(double *lbox_1, double *lbox_2, int npart_1, int npart_2, 
		  char *file_1 , char *file_2, double dz){
  int dummy;
  double r[3], v[3], m, z, lbox_new[3];
  char type;

  for ( int k=0; k<2; k++ ){
    if ( lbox_1[k] < lbox_2[k] ) lbox_new[k] = lbox_2[k];
    else lbox_new[k] = lbox_1[k];
  }
  lbox_new[2] = lbox_1[2] + lbox_2[2] + dz;


  FILE *fin1 = fopen(file_1, "r");
  FILE *fin2 = fopen(file_2, "r");
  FILE *fout = fopen("start.xyz", "w");

  if ( fin1 == NULL || fin2 == NULL || fout == NULL )
    merror("%s at %d: Couldn't open files", __func__, __LINE__);

  fprintf(fout, "%d\n", npart_1+npart_2);
  fprintf(fout, "%f %f %f\n", lbox_new[0], lbox_new[1], lbox_new[2]);
  

  // copy first file contents
  fscanf(fin1, "%d\n", &dummy);
  fscanf(fin1, "%lf%lf%lf\n", &r[0],&r[1], &r[2]);
  
  for (int n=0; n<npart_1; n++ ){
    fscanf(fin1, "%c%lf%lf%lf%lf%lf%lf%lf%lf\n", &type, &r[0],&r[1],&r[2], 
	   &v[0],&v[1],&v[2], &m,  &z);
    fprintf(fout, "%c %f %f %f %f %f %f %f %f\n", type, r[0],r[1], r[2], 
	    v[0],v[1],v[2], m,  z);
  }

  fclose(fin1);

  // second file with shift
  fscanf(fin2, "%d\n", &dummy);
  fscanf(fin2, "%lf%lf%lf\n", &r[0],&r[1], &r[2]);
  
  double shift = dz+lbox_1[2];
  for (int n=0; n<npart_2; n++ ){
    fscanf(fin2, "%c%lf%lf%lf%lf%lf%lf%lf%lf\n", &type, &r[0],&r[1],&r[2], 
	   &v[0],&v[1],&v[2], &m,  &z);
    fprintf(fout, "%c %f %f %f %f %f %f %f %f\n", type, r[0],r[1], r[2]+shift, 
	    v[0],v[1],v[2], m,  z);
  }

  fclose(fin2);

  fclose(fout);

  printf("Wrote xyz configurational file\n");
}
